Mathematical reconstruction of the metabolic network in an in-vitro multiple myeloma model

It is increasingly apparent that cancer cells, in addition to remodelling their metabolism to survive and proliferate, adapt and manipulate the metabolism of other cells. This property may be a telling sign that pre-clinical tumour metabolism studies exclusively utilising in-vitro mono-culture models could prove to be limited for uncovering novel metabolic targets able to translate into clinical therapies. Although this is increasingly recognised, and work towards addressing the issue is becoming routinary much remains poorly understood. For instance, knowledge regarding the biochemical mechanisms through which cancer cells manipulate non-cancerous cell metabolism, and the subsequent impact on their survival and proliferation remains limited. Additionally, the variations in these processes across different cancer types and progression stages, and their implications for therapy, also remain largely unexplored. This study employs an interdisciplinary approach that leverages the predictive power of mathematical modelling to enrich experimental findings. We develop a functional multicellular in-silico model that facilitates the qualitative and quantitative analysis of the metabolic network spawned by an in-vitro co-culture model of bone marrow mesenchymal stem- and myeloma cell lines. To procure this model, we devised a bespoke human genome constraint-based reconstruction workflow that combines aspects from the legacy mCADRE & Metabotools algorithms, the novel redHuman algorithm, along with 13C-metabolic flux analysis. Our workflow transforms the latest human metabolic network matrix (Recon3D) into two cell-specific models coupled with a metabolic network spanning a shared growth medium. When cross-validating our in-silico model against the in-vitro model, we found that the in-silico model successfully reproduces vital metabolic behaviours of its in-vitro counterpart; results include cell growth predictions, respiration rates, as well as support for observations which suggest cross-shuttling of redox-active metabolites between cells.


Introduction
A cell's metabolism can be defined by extensive interconnected chemical reaction networks subdivided by seemingly organised pathways [1,2].These pathways specialise in building unique blocks necessary for all aspects of cellular functioning (i.e., homeostasis, cell reproduction, and biomass) [3,4].Their performance depends on many aspects, a critical one being the microenvironment in which the cell resides [4][5][6][7].A microenvironment consists of a reservoir of readily available chemicals and different cell types.As the microenvironment's cell population have different metabolic requirements, their spatial organisation, type, and objectives will be affected by the properties of the microenvironment itself.This is particularly evident in the case of cancer metabolism [8].Studies have demonstrated that several types of malignancies rely on their immediate extracellular environment and the cells that reside within to meet their highly energetic demands, enhance proliferation, and even enable therapy resistance [5][6][7][9][10][11][12][13].
A powerful tool for probing the metabolic interactions between cell types are in-vitro cocultures (Fig 1) [10][11][12][13].These experimental models are advantageous as cell interactions can be studied under a controlled microenvironment that allows the co-cultured cells to establish a unique metabolic niche [14][15][16].The resulting intercellular chemical exchanges enable cells in the co-culture to tap into a larger pool of metabolites and, thus, instigate a substantial division of metabolic labour that aims to optimise the use of available nutrients.Consequently, cell survival and proliferation may be enhanced, particularly under system perturbations (e.g., hypoxia or lack of oxygen) in which mono-cultures would otherwise perish [14][15][16].
Co-culture studies are possible today thanks to the technological advancement seen during the past couple of decades [17].Nevertheless, the study of intercellular metabolic interactions across microenvironments remains incredibly challenging [18][19][20].As a result, a wave of interdisciplinary research avenues has given rise to numerous bioinformatic approaches that unravel cell metabolism.Techniques like omic enrichment and extensive multi-scale and multi-omics analyses have helped underpin critical aspects of cell functionality with unprecedented precision [21][22][23].However, the granularity required to fully understand the complexities of cell biochemistry, such as quantitative information on metabolic fluxes, is, more than often, too challenging to obtain, analyse, and interpret.This is an all-too-familiar problem in multiple myeloma (MM) research.MM is a type of malignancy found in the bone marrow (BM) and is formed through the ingress and proliferation of malignant plasma cells [11,24].The disease remains incurable, and even though novel therapies have improved patient survival outcomes, more effective treatments are needed for many patients unable to tolerate current cancer management strategies [25,26].Therefore, developing next-generation therapies for multiple myeloma requires a technical approach focusing on delivering long-term disease control while maintaining the quality of life [26].
A potential approach to address this issue involves exploiting those metabolic features critical for MM survival.Malignant plasma cells have been observed to transform the phenotype of the stroma residing in the bone marrow niche [7,15,27].This malignancy can generate a metabolic network within the tumour microenvironment capable of promoting survival, and aggressive proliferation [28,29].Studies have previously demonstrated that the metabolism of the bone marrow is significantly altered in patients with MM and that the bone marrow mesenchymal stem cell (BMMSC) is a primary supportive cell type for malignant plasma cells [28,30].We thus hypothesise that understanding the physiology of the metabolic interaction between these two cell types (the BMMSC and the MM cell) and the niche in which they reside will facilitate the discovery of a novel metabolic target.Such a target could interfere with a metabolic intercellular pro-survival axis presumed to be critical for the malignancy's survival.
In this light, our study outlines a multidisciplinary approach to produce a testable and integrated in-vitro/in-silico model of the metabolic network formed by malignant plasma cells, or myeloma, and bone marrow mesenchymal stem cells.We achieve this via two phases that display the synergy between computational and bench-side experimental research.The first or experimental phase involves mono-culturing HS-5 (BMMSC) and JJN-3 (MM) cell lines.We performed benchmark measurements such as cell growth and cell respiration rates, that is, oxygen (O 2 ) consumption, to probe their carbon metabolism capacity [31,32].We followed these experiments by in-vitro co-culturing the aforementioned cell lines (Fig 1A).Finally, we probed the metabolic phenotype of the co-culture system using stable-isotope 13 C 6 Glucose and 13 C 5 Glutamine tracing [33,34].
The computational phase of this study draws from and optimises established techniques to generate two cell-specific genome-scale models (GEMs) of our two cell types, the BMMSCs, and the MM [35,36].Using these GEMs together, we reconstruct an in-silico co-culture constraint-based model (CBM).The model is assembled through a bespoke computational workflow that integrates publicly available transcriptomic data sets (retrieved from the NCBI Gene Expression Omnibus-GEO), using the mCADRE algorithm, to a thermodynamically constrained version of the latest global knowledge base of metabolic functions categorised for the human genome: the Human Recon3D; using a subroutine of the redHuman algorithm (Fig 2) [37][38][39][40].Our computational co-culture was further constrained to mimic the in-vitro growth medium, RPMI, using the metabotools MATLAB routine [41].We validate the resulting model against its experimental counterpart, where we found it successfully recapitulated the observed in-vitro growth and respiration rates.Furthermore, via integration of stable-isotope labelling data via 13 C-metabolic flux analysis ( 13 C-MFA), using the MATLAB INCA routine as a means of "phenotype tuning" [42,43].

Cell culture
The HS-5 cell line (ATCC) and JJN-3 cell lines were both STR profiled to ensure accurate identification and maintained in RPMI (Sigma-Aldrich, R8758) with 10% FBS (Sigma-Aldrich, F7524).For individual cell growth experiments, HS-5 and JJN-3 cells were seeded at 4 × 10 4 cells/mL, and 3.7 × 10 4 cells/mL and cell growth was assessed over time by cell count for JJN-3 and by sulforhodamine B (SRB) assay for HS-5.For co-culture experiments, HS-5 cells were seeded at 12 × 10 4 cells per well of a 6-well plate the day before JJN-3 were seeded on top of a transwell (CC401, Appleton Woods) at 8 × 10 4 cells per well.

Thermodynamic data integration
Gibb's Free Energy De nes Directionality

Oxygen consumption
Trypsinized cells were resuspended in culture media and loaded in an Oxygraph-2k (Oroboros instruments) chamber.After closing the chambers and recording routine respiration, Oligomycin (2.5 μM) was added to inhibit ATP synthase.Respiration was inhibited by the addition of Rotenone (0.5 μM) and Antimycin A (2.5 μM) at the end of the experiment.Measurements of the non-phosphorylating electron transfer system (ETS) capacity were obtained through stepwise (0.5 μM) titration of the uncoupler, carbonyl cyanide 4-(trifluoromethoxy) phenylhydrazone (FCCP).

Carbon metabolic tracing
Cells at 70% confluency were incubated in flux media: modified RPMI 1640 without glucose, glutamine, L-isoleucine, L-leucine, L-valine, and phenol red (Cell Culture Technologies) supplemented with either 2g/L 13 C-[U]-glucose or 2mM 13 C-[U]-glutamine (both CK isotopes), with the other provided in an unlabelled form.After 48 hours, media was removed for extraction, JJN-3 cells spun down, and then cells (JJN-3) and wells (HS-5) were washed twice with ice-cold saline before quenching metabolism with ice-cold MeOH.Cells were transferred to a cold tube into which was added D 6 -glutaric acid in ice-cold H 2 O (CDN isotopes, D-5227) and pre-chilled chloroform.After shaking on ice for 15 minutes and centrifugation, the polar phase was transferred to another tube to be dried.

Derivatisation and GC-MS
Dried-down extracts were derivatized with 2% methoxamine in pyridine (20μl, for one hour at 60˚C) followed by N-(tert-butyldimethylsilyl)-N-methyl-trifluoroacetamide, with 1% tert -butyldimethylchlorosilane (30μl, 1h at 60˚C).Samples were transferred to glass vials for GC-MS analysis using an Agilent 8890 GC and 5977B MSD system.1μL of the sample was injected in splitless mode with helium carrier gas at a 1.0 mL/min rate.The initial GC oven temperature was held at 100˚C for one minute before ramping to 160˚C at a rate of 10˚C/min, followed by a ramp to 200˚C at a rate of 5˚C/min and a final ramp to 320˚C at a rate of 10˚C/ min with a five-minute hold.Compound detection was carried out in scan mode.Total ion counts of each metabolite were normalised to the internal standard D 6 -Glutaric acid.

Normalisation and quantification
Data normalisation occurred by counting cells using counting chambers Fast Read 102, purchased from Kova International.GCMS data were analysed using Agilent Mass Hunter software for real-time data quality analyses and in-house MATLAB scripts.

Model
A frequently sought-after solution for producing in-silico large metabolic models is to map out tissue-specific functions using genome-scale models (GEMs) [44,45].GEMs are comprehensive models that aim to recreate a particular organism's metabolism [36,38,44,45].They typically include most, if not all, known chemical reactions and metabolites in the form of a stoichiometric matrix and their corresponding genes in the form of a gene-rule matrix.The latter is an attractive feature of GEMs because each enzyme-mediated reaction in the model is associated with information regarding how a gene influences a protein, which in turn influences a reaction.In GEMs, this is referred to as a gene-protein-reaction (GPR) rule.GPRs are essential to GEMs as their description of the relationship between genes encoding enzymes enables the mediation of a given reaction or set of reactions in a given model [44,46].
Furthermore, GEMs are attractive because they promise to advance our understanding of the underlying metabolism behind various physiological and pathological processes in various tissues [44,46,47].Unlike many modelling frameworks, GEMs can yield a profound insight regarding cell behaviour requiring minimal information on the biophysical equations that require difficult-to-measure kinetic parameters.Efforts such as the global knowledgebase of metabolic functions categorised for the human genome, known as Human Recon3D, coupled with abundant high-throughput data, make the reconstruction of cell-specific metabolic models possible [39,48].What follows is a detailed workflow of a top-down reconstruction strategy designed to produce two cell-specific constraint-based models; one bone marrow mesenchymal stem cell (BMMSC) and one myeloma (malignant plasma-MM) cell (Fig 2).Our overarching aim is to assess the consequences of metabolic adaptation and regulation of the desired phenotype-in our case, a multiple myeloma in-vitro co-culture; via network flux analysis.To this end, the first step uses the generic model of human metabolism: Recon3D, the latest and arguably the most comprehensive consensus human GEM.Recon3D consists of 10,600 reactions, 5938 associated with 2248 genes and 2797 unique metabolites across seven compartments: cytosol, mitochondria, peroxisome, Golgi apparatus, endoplasmic reticulum, nucleus, and lysosome [39,48].GEMs of metabolic networks combined with CBM produce a powerful tool that ultimately yields qualitative and quantitative information regarding the distribution of steadystate metabolic fluxes through a given network graph representation.This feature enables the prediction of experimentally measurable variables such as cell growth rates, or production rates of metabolites of interest [49].In this paradigm, the steady-state approach works because, within the time scale of our experiments, changes in the concentration of metabolites occur slowly [49][50][51].This is even true under constant exponential cell growth (i.e., cancer).The steady-state assumption allows CBMs to model large networks without knowing enzyme kinetics and post-translational regulatory mechanisms [50].This last feature is advantageous when cellular activity in interest physiology is not well-understood [50,52].

Thermodynamic data integration (Step 1)
We introduced thermodynamic information for all metabolite compounds and reactions by following the protocol outlined in the redHuman workflow [40].Briefly, the protocol involves obtaining metabolite Gibb's free energy data from the MetaNetX archive [53].Those metabolites from Recon3D with identifiers from SEED, KEGG, CHEBI, and HMDB, are manually annotated and then, using ChemAxon's Marvin (available for free if used for academic purposes), the compound structures were transformed into their primary protonation states at a pH of 7 and generated MDL Molfiles [53][54][55][56][57][58][59].These MDL molfiles and the Group Contribution Method were then used to estimate the standard Gibb's free energy of the formation of the metabolites in Recon3D as per [60,61].As a result, each metabolite is associated with its Gibb's free energy (ΔG).Then the ΔG of a given reaction is estimated from the thermodynamic properties of its reactants and products.Following the incorporation of this data, the model was then primed for thermodynamic flux balance analysis (tFBA) and thermodynamic flux variability analysis (tFVA) using the prepModelforTFA.mroutine from redHuman [40].A readily available annotated matrix and a MATLAB script are available in the matTFA package by [62].We, however, repeated the protocol due to updates in the MetaNetX archive (with the latest update occurring in March 2022).The model was verified for consistency by running the redHuman's matTFA subroutine and comparing its results to the thermodynamically annotated Recon3D version generated by the redHuman algorithm package [40,62].
Performing this step allows one to determine and ensure the feasible thermodynamical directionality of those chemical reactions included in the Recon3D matrix for which we have no constraining data.Numerically, this also ensures that the employment of optimisation routines, such as optimiseCBmodel from the MATLAB COBRA routine, does not produce thermodynamically infeasible flux loops, thus avoiding physiologically inaccurate behaviour [63].

Database mining and cell-specificity (Step 2)
To determine which reactions must be removed from the now thermodynamically constrained Recon3D matrix to generate our cell-specific models, we used transcriptomic datasets from the Gene Expression Omnibus database (GEO) [37,39].In this study, we strictly adhere to those data sets based on the Affymetrix Human Genome U133 Plus 2.0 Array platform for the bone marrow mesenchymal stem cell and the malignant plasma cell models [64].The choice was driven by full support from Affymetrix via the availability of the .cdffile (a file that describes the layout for an Affymetrix GeneChip array) and its compatibility with the latest MATLAB version through the Bioinformatics Toolbox [64].Those sets for generating the bone marrow mesenchymal stem cell model reconstruction were obtained from healthy human early bone marrow passage samples and gene expression data from mesenchymal stem cells cultured in MSCGM (mesenchymal stem cell growth medium).These were chosen in line with an early model that uses a legacy version of Recon3D, namely Recon1, to generate a genome-scale reconstruction-iMSC1255 [65,66].Their accession ID: GSM184636, GSM184637, GSM184638, GSM194076, GSM194077, GSM194078, GSM194079, GSM764199, GSM797497, GSM797498, GSM920586, GSM920587, and GSE80608 [67].In particular, the GSE80608 data set includes newer datasets gather from gene expression of bone marrow mesenchymal stem cells grown out from bone marrow aspirates for multiple patient passages [30].
For the myeloma cell model, we have used those datasets from a series of pre-treatment bone marrow aspirates derived from patients diagnosed with multiple myeloma.The myeloma cells used in these samples were selected for their expression of CD138 + , a characteristic marker of multiple myeloma.The data sets also contain expression of individuals presenting a condition called Monoclonal Gammopathy of Unknown Significance (MGUS).Athough this condition has been found to represent a precursor to multiple myeloma, its important to note the significant clinical difference between these patient groups.Hence, in this study, we did not use data from individuals with MGUS, acknowledging the low probability transition rate from MGUS to multiple myeloma (1% annually) [30].The data accession IDs: GSM50986, GSM50987, GSM50988, GSM50989, GSM50990, GSM50991, GSM50992, GSM50993, GSM50994, GSM50995, GSM50996, and GSM50997 [68].It is important to remark that the expression data above provides a rich source of information about the cellular environment.However the model tuned with this data alone may not fully replicate in-vivo or in-vitro microenvironments.Consequently, the inclusion of this gene expression data may cause alterations when compared to experimental results.In other words, if the model were to be informed solely by early human bone marrow samples, differences might arise due to the specific characteristics of these samples, which can vary from the conditions in growth media.The inclusion of different data types is necessary for our study (see "Phenotype specification" section below).
Using the thermodynamically constrained Recon3D (Step 1) and the microarray datasets above, we executed the mCADRE routine [38].Briefly, as a first step, mCADRE scores those genes included in the Recon3D matrix according to their presence in the transcriptomic dataset.These scores are then attributed to reactions according to the initial model's annotated gene-protein-reaction (GPR) relations.Next, the topology of the metabolic network is considered to update these scores.Once the procedure is complete, reactions with scores below a certain threshold are removed to complete the reconstruction of the cell-specific metabolic network.We refer our reader to [38] for a detailed explanation of the mCADRE algorithm.
We have devised an updated version of mCADRE, compatible with Recon3D (https:// github.com/esig626/Metabolic_Networks).This version operates with a Matlab wrapper to enable the use of IBM's CPLEXstudio Linear Programming (LP) and MLP solvers directly in their native C version.This is necessitated due to IBM's discontinuation of support for the MATLAB software in the most recent verisons [38,69].Moreover, despite not being open source, IBM's CPLEXstudio can be accessed for research purposes under academic and student licenses.This accessibility is vital, given the mCADRE algorithm's reliance on CPLEX's fast flux variability analysis, also known as fastFVA, optimisation routine [70].In light of these considerations, we have incorporated an optimised, open-source version of fastFVA, aptly named Very Fast Flux Variability Analysis (VFFVA).VFFVA leverages CPU parallelisation, thereby facilitating biological insight generation through optimised system biology modelling.It is available, under an open source license, for use in C, MATLAB, and Python [71].

Model reduction (Step 3)
In the previous steps, we have outlined how our GEMs, the BMMSC and Myeloma, were reconstructed to study the biochemistry occurring in human cells.In essence, these models can be used to probe the biochemistry of their respective metabolism.However, the validation of this approach's results is limited because most experimental procedures focus on a small subset of cell metabolism.As a result, the model and its solutions leverage unnecessary complexity in its analysis, hindering a consistent and concise physiological representation.In our study, we are interested in the metabolism of cancer cells, and the data we can produce is centred around central carbon metabolism; this means we must keep those pathways that provide the energy, redox potential, and biomass precursors for cell growth and sustenance.We also include those pathways reported to be altered in cancer cells.Consequently, we specify the core subsystems of our models to be glycolysis, the pentose shunt, tricarboxylic acid cycle, oxidative phosphorylation, glutamate/glutamine metabolism, serine metabolism, urea cycle, and reactive oxygen species detoxification.We executed the redHuman algorithm, which performed a series of reductions [40].Briefly, the core metabolic pathways defined above were specified in the routine.We then used the mCADRE GEM as input, which redHuman uses its subroutines, redGEM-lumpGEM, to reconstruct a reduced model around the core we specified [38,40,72,73].

Co-culture assembly, objective function, medium and phenotype specification (Step 4)
Co-culture assembly and objective function.To model the experimental cell culture conditions, we must define an in-silico culture medium based on those used in the in-vitro counterpart.A correct definition of the medium in the model is fundamental for adequately representing intracellular metabolism.Using a combination of the redGEMX (part of redHuman) and the metabotools MATLAB package routine setMediumConstraints, we could guarantee that the reduced model has all the feasible pathways that consume and produce the components of the extracellular environment of the medium [40,41,72].The redGEMX subroutine finds those pathways from the model needed to connect the extracellular metabolites to the core network defined above.We also employed constraints on the exchange fluxes to the extracellular medium in both cells by using data related to RPMI media composition.The advantage of using these methods in conjunction is that we were able to specifically tend to the needs of the reduced model whilst providing a physiologically sound in-vitro medium model that defines the medium components and those that can be uptaken by the cells but are not captured by the measured data such as CO 2 , H + , H 2 O, HCO À 3 , and NH 4 .When considering the assembly of the in-silico co-culture model, we manually curated each cell-specific stochiometric matrix to reflect the compartments of each as unique (i.e., metabolite[cm] for myeloma cytosol and metabolite[bm] for BMMSC cytosol), excluding the extracellular fluxes which remain the same given that both cells share this space (specified in the model under the notation [e]).The consequent model spans a large metabolic matrix where duplicated extracellular sinks and sources were removed.Furthermore, the final model requires that both cells are equipped with biomass-producing reactions.We adhere to those specified by the Recon3D GEM [39].There are several reasons, including experimental capabilities, for using this function as the optimisation objective of the model.One such argument was brought forward by a study which analysed and compared several CBM biomass functions for cancer metabolic modelling [4,74].Here, it was found that the formulation from the Recon family models (Recon 2 and Recon3D) leads to the most accurate results in terms of essential gene and growth rate predictions amongst all cancer cell lines they tested, including those for leukaemia research.
Multi-objective optimisation.Our two cells residing within the microenvironment (i.e., the extracellular media) have different functions, each needing a different objective in a simulation.In a multi-objective paradigm, we assume that both cells, even under controlled environments, address various biological tasks (i.e., objectives) to sustain life.To obtain an optimal solution, we use the concept of "Pareto-optimality".Briefly, a Pareto-optimal solution is one where the performance of one task would diminish the ability to achieve one or more other tasks.As we assume, for simplicity's sake, that regardless of the biological functions of the two cell types, in our simulation, the biomass function carries sufficient knowledge about the biological system [75][76][77].Our simulations implemented multi-objective optimisation using the MOFA routine for the COBRA package in MATLAB [78].Briefly, MOFA maps the n-1-dimensional surface of the Pareto front, where n is the number of objectives, by calculating a large set of Pareto solutions.That is, the routine maps an area of multiple criteria decisionmaking in the presence of trade-offs between two or more conflicting objectives [78,79] (S1 Text).
Upon curation of the co-culture model, we performed flux variability analysis (FVA), and used the average between the extremes of minimisation and maximisation of the objective function.an extension of Flux Balance Analysis (FBA).Briefly, FVA explores the range of possible flux values for each reaction in a metabolic network, given a certain objective function.In doing so, we not only ensured a robust flux distribution but also fostered a superior prediction, all while preserving the essential functionality of pathways integral to biomass synthesis [71,80,81].However, it must be noted that the effectiveness and accuracy of this approach depends on several factors, including the quality and completeness of the metabolic model, the choice of objective function, and the accuracy of the input data.The method also provides an exhaustive window of exploration of the cellular metabolic network.It allows us to mitigate potential bottlenecks that could arise from overconstrained solutions, such as the over-maximisation or minimisation of the objective function.Concomitantly, it enables the identification of potential metabolic fluxes of interest.In our model, the inherent complexities of the metabolic reactions that remain consistent despite flux variations and inherent uncertainties often function suboptimally due to various constraints and trade-offs.This is particularly notable in complex biological systems like mammalian cells, which demand a fine balance between growth rate maximisation, energy conservation, and redox balance maintenance as objectives [71,81].
By providing a more nuanced representation of the complex interplay between diverse metabolic pathways and their interactions, screening the metabolic network in this manner offers a comprehensive understanding of the processes that underpin the co-culture system.(S3 Table ).
Phenotype specification.Because metabolism is highly cell type-and state-specific, it offers valuable opportunities for diagnosing and treating cancers [82].Mass spectrometry (MS) or nucleic magnetic resonance (NMR) based methods such as metabolomics measure relative abundance of cellular metabolites, which can help detect differences in metabolic state, but such abundance data does not inform on itself metabolic activities [42,83,84].For a detailed assessment of intracellular metabolic activity in cells, stable isotope labelled nutrients are required [83][84][85].This is because metabolic reaction rates, or fluxes, contribute to metabolic phenotypes and mechanisms of cellular regulation [42,83,85].However, in most cases, isotopic labelling data cannot be interpreted intuitively due to the highly complex nature of atom rearrangements in metabolic pathways; instead, a formal model-based analysis approach is required to extract flux information from the labelling data.To quantify and characterise the metabolic phenotype in cells 13 C-metabolic flux analysis ( 13 C-MFA) is the gold standard for converting isotopic labelling data into corresponding metabolic flux maps [42,43].At a highlevel 13 C-MFA is formulated as a least-squares parameter estimation problem, where fluxes are unknown parameters that must be estimated by minimising the difference between the measured labelling data and labelling patterns simulated by the model, subject to stoichiometric constraints resulting from mass balances for intracellular metabolites and metabolite labelling states; the so-called mass isotopomers [42,43].
Our main objective is to generate a quantitative map of cellular metabolism by assigning flux values to the reactions in an annotated metabolic network model and confidence intervals for each estimated flux.The confidence intervals act as lower and upper bounds in the CBM generated from steps 1 through 4. The result is the incorporation of experimental data and further modelling constraints.Consequently, we allow each cell in the model to acquire the phenotype displayed, or extrapolated, from the co-culture in-vitro setup.Although many routines can perform 13 C-MFA, we have chosen the MATLAB package for isotopomer network compartmental analysis (INCA) [43].Briefly, INCA uses the elementary metabolite unit (EMU) framework to simulate isotopic labelling in any arbitrary biochemical model.The routine is equipped with a Montecarlo parameter estimation routine, which we used to obtain confidence intervals for each predicted flux [42,43,86].
The selection of isotopic tracers in this study was 13 C 6 -Glucose and 13 C 5 -Glutamine. 13C 6 -Glucose is best for determining those fluxes in the upper sections of carbon metabolism, that is, glycolysis and the incorporation of its products into the TCA cycle [42]. 13C 5 -Glutamine, on the other hand, is used to determine exclusively lower parts, that is, the TCA cycle and reductive carboxylation.All tracing experiments performed were done in parallel, meaning that the cultures and initial conditions were identical, except that tracing was done in different batches.The data obtained from these experiments can also be incorporated in parallel because these isotopes are complementary [42,43] (S2 Text and S2 to S8 Tables).

Models in mono-culture: Results and experimental validation
We performed benchmark tests to validate the predictive behaviour of both models, the BMMS cell and the MM cell.The following simulations were performed in mono-culture.
Each cell-specific model was tested individually, where the extracellular medium was set to RPMI.We then performed, in-silico, experiments that address the consistency of cellular cultures [87,88] (S1 Table ).
Cell growth.Simulation for cell growth prediction rate was obtained by performing the Flux Balance Analysis routine with maximisation of biomass as the optimisation objective [48].BMMSC growth rate has been reported to largely depend on cell culture conditions, including medium composition and its changes over time.Our experimental growth rate data from the HS-5 cell line cultured in RPMI medium was calculated to be 0.024 ± 0.0097 /hour.Several works of literature report the in-vitro BMMSc population growth rate to be 0.027 /hour.Our bone marrow mesenchymal stem cell model predicts a growth rate of 0.02367 /hour.This was found by noting that each gram of dry cell weight (grDW) is assumed to be equivalent to one millimole of biomass.This means that the optimal in-silico specific growth rate is equal to [65]: Our results are within the measured error of our experiments and those reported in literature (Fig 3A) [65,[89][90][91].We repeated the procedure to calculate the growth rate of the myeloma cell (Eq 1).This growth curve is contrasted with the projected growth as determined by the exponential growth model that uses the growth rate predicted from the BMMSC Constraint-Based Model (CBM); solved utilising the Flux Balance Analysis (FBA) procedure with the maximum biomass production rate serving as the optimization objective.This is represented by the black intermittent line.Simultaneously, the figure demonstrates the measured, normalised growth of the MM cell line in-vitro culture, symbolized by red dots.This is set against the growth curve predicted, with maximum biomass production rate being the primary optimisation goal; denoted by the intermittent red line.B. Illustration of the growth rate comparison between the model and experimental results in co-culture.The figure demonstrates the measured, normalised in-vitro co-culture growth of the BMMSC cell line; (black dots).This is contrasted with the growth curve given by the exponential growth model solution, which is fitted with the growth rate predicted by the in-silico co-culture model, represented by the intermittent black line.Concurrently, the figure presents the measured, normalised in-vitro co-culture growth of the MM cell line, depicted by red dots.This is set against the growth curve yielded by the exponential growth model solution, adjusted to the growth rate forecasted by the in-silico co-culture model, denoted by the intermittent red line.Note that in both cases, the monoculture and the co-culture, the negligible differences between the modelled and measured growth rates of BMMSCs and MMs limit cell type discrimination. https://doi.org/10.1371/journal.pcbi.1011374.g003

Oroboros respiration study.
Because carbon metabolism is central to cancer study, as described in the methodology section, we tested each model's (i.e., mono-culture) respiration by simulating an in-silico Oroboros study [32].We then compared those performed in the invitro culture (also in mono-culture).Briefly, an Oroboros respiration study consists of a series of electron transport chain perturbations via pharmacological inhibitors to measure the respiratory capacity of the cell [32].The procedure commences after basal measurements are recorded.This was achieved by measuring the oxygen (O 2 ) flux when performing an FBA optimisation routine, where the objective function (Biomass) was set to find the minimum flux distribution through the network (Fig 4A and 4B α, purple and green; respectively; Fig 4A and  4B α, red and yellow; respectively).The second step in the respiration study involves exposing cells to Oligomycin, an inhibitor of the ATP synthase (or complex V of the electron transport chain).By inhibiting the proton (H + ) flux through this enzyme, the effects of Oligomycin cause an increased proton gradient across the mitochondrial inner membrane preventing electron transport through complexes I-IV (Fig 4A and 4B β purple and green; respectively, and Fig 4A and 4B β, red and yellow; respectively).Oxygen consumption falls proportionally.The remaining rate of mitochondrial respiration represents a proton leak; that is, the protons pumped during electron transport that result in oxygen consumption but not ATP production (Fig 4A and 4B.δ, purple and green; respectively, and Fig 4A and 4B δ red and yellow; respectively).An increase in ATP-linked respiration, a measure of the cell's capacity to meet its energy demands, would indicate an increase in ATP demand.In contrast, a decrease would indicate either low ATP demand, a lack of substrate availability, including severe damage to the oxidative phosphorylation pathway, which would impede the flow of electrons and result  4B �, red and yellow; respectively).In both models, we achieved this by constraining the ATP synthase flux to zero and constraining the electron transport chain complexes to prevent reversal.Then we performed a second FBA optimisation routine, where the objective function was to maximise the biomass production rate.
The following experiment involves FCCP, a potent uncoupling agent that collapses the proton gradient and disrupts the mitochondrial membrane potential.This agent is added following the cell's exposure to Oligomycin.As a result, electron flow through the electron transport chain is uninhibited, and oxygen consumption by complex IV reaches the maximum.A high FCCP-stimulated oxygen consumption rate compared to a basal rate indicates that the mitochondria are using less than the maximal rate of electron transport supported by substrate supply from the cells.We performed this simulation by relaxing the constraints in the electron transport chain-related fluxes.We then performed a third FBA optimisation routine, where the objective was to maximise the flux through complex IV of the model's electron transport chain (Fig 4A and 4B γ, purple and blue; respectively, and Fig 4A and 4B γ, red and yellow; respectively).We then used the FCCP-stimulated oxygen consumption rate to calculate the cell's spare respiratory capacity.This rate is defined as the difference between maximal respiration and basal respiration.Spare respiratory capacity measures the cell's ability to respond to increased energy demand or under stress.
The last experiment involves exposing cells to a mixture of Rotenone, a complex I inhibitor, followed by the addition of antimycin A, an electron transport chain complex III inhibitor.This combination shuts down mitochondrial respiration and enables the calculation of nonmitochondrial respiration driven by processes outside the mitochondria.We performed this insilico experiment by constraining the flux through complexes 1 and 3 of the electron transport chain.We performed an FBA optimisation routine where biomass was used as the objective function to maximise (Fig 4A and 4B z, purple and green; respectively, and Fig 4A and 4B z, red and yellow; respectively).

Models in co-culture: Results and experimental validation
Co-culture cell growth.Following the experimental validation of the mono-culture models, we verified the co-culture in-silico model growth rate prediction.We found its value to be 0.0335 /hour for both cells, the BMMSCs and MM cells alike (Fig 3B).Interestingly, although it appears remarkably close to the measured values in its in-vitro counterpart, we cannot accurately validate the growth rate for the JJN-3 cell line in co-culture.This is because our experimental data was only performed for one biological replicate, rendering its statistical significance null.Because the experimental data appears inconclusive, we assume the model's prediction suffices as a reasonable approximate growth according to those found in mono-culture, and those seen in literature [65,[89][90][91][92][93][94].
13 C-MFA.Mass and isotopomer balances were simulated for Myeloma cell and Bone Marrow Stem Cell, HS-5 and JJN-3 cell lines.Isotopic steady state was verified by measuring mass isotopomer distributions (MIDs) equilibration at the time points analysed.The isotopic steady state is reached when the change in isotopic enrichment over time falls within the measurement uncertainty range.We implemented the 13 C 6 Glucose and 13 C 5 Glutamine datasets in parallel datasets by regressing simultaneously to yield one complete metabolic flux map for each medium formulation.Glucose uptake flux was calculated to be 0.45 μmol/10 6 Cells/hr, and the measured glutamine was 0.63 μmol/10 6 Cells/hr-these values were calculated according to the protocol outlined in [42].All fluxes were simulated in the INCA routine for MATLAB using a minimum of 100 unique restarts from random initial values to ensure a global minimum was found [43].Flux results were subjected to a chi-square statistical test to assess goodness-of-fit, and 95% confidence intervals were calculated for each estimated flux value [43].illustrate the observed and simulated mass isotopomer distributions (MIDs); also known as "labelling patterns."When glucose is metabolised, the incorporation of isotopes leads to a shift in mass (m+n, where n ranges from 1 to 5 in Fig 5).These MIDs represent the relative abundance of each isotope, normalised to the total sum of all possible isotopologues.Hence, a  13 C]Glucose and [U- 13 C]Glutamine tracing.The Mass Isotopomer Distributions represent the relative abundance of each isotope, normalised to the total sum of all possible isotopologues.Glucose-derived metabolite abundances.Data suggests high glycolytic activity in both cells.In JJN-3 cells, pyruvate-derived alanine is much higher than in HS-5 cells.The abundance of pyruvate in JJN-3 cells is also more significant than that of HS-5.It is possible to see that although pyruvate-derived citrate follows the expected TCA cycle, it does not continue downstream as aketoglutarate is low in JJN-3 cells.This indicates an interruption of the TCA cycle.C. & D. Mass isotopomer distributions from [U- 13 C]Glucose and [U- 13 C] Glutamine tracing.Glutamine-derived metabolite abundances.Data suggests that the TCA cycle in both cells is supported via anaplerosis, where exogenous glutamine is interconverted to a-ketoglutarate and then incorporated into the TCA cycle.https://doi.org/10.1371/journal.pcbi.1011374.g005metabolite containing n carbon atoms may have 0 to n of its atoms labeled with 13 C, resulting in isotopologues that exhibit an increase in mass from m+0 to m+n for a specific metabolite.For an in-depth explanation of the concepts used in 13 C-MFA we refer the reader to antoniewicz et.al, (2018) [42].These results are thoroughly discussed in the next section (S2 Text and S2 to S7 Tables).
Co-culture model predictions.Our in-silico simulation results reveal that both cells decrease their glycolytic activity as they go from mono-culture to co-culture.This, in principle, is a counterintuitive result as a high glycolytic activity is a hallmark of cancer metabolism.Therefore, we expected to see an enhancement of this behaviour when assembled in a more physiologically relevant way.However, this result capitulates the observed heterogeneity between cell types in mono-culture, as opposed to its results in co-culture.For instance, the mono-culture solutions of our models prioritise those fluxes associated with the glycolysis pathways to generate most of the ATP required for biomass synthesis.In contrast, we see this to a lesser extent in in-vitro co-cultured cell-lines.Contrary to our expectations, our mono-culture models support this experimental evidence where the in-vitro co-culture cell lines consume less glucose than mono-cultured.However, glucose uptake simulations performed in our co-culture model do not re-capitulate this behaviour, perhaps explaining the slightly faster cell growth rate result observed in Fig 3B .The glucose usage phenotype that the in-vitro co-culture exhibits was acquired once the 13 C-isotope labelling data was incorporated.It has been suspected that synthesising large amounts of bioenergetic metabolites requires malignant cells to establish a cooperative intercellular metabolic network.Our model predicts that a critical component of this mechanism is driven by the unidirectional export of pyruvate and lactate from bone marrow mesenchymal stem cell to the myeloma cell.Our model predicts that this now exogenous pyruvate might be uptaken by the myeloma cell to support its respiratory and biomass needs.In Fig 6, it is depicted how in co-culture, BMMSC produces relatively large amounts of pyruvate relative to those it uses for its functionality.Indeed, a fraction of this branching point metabolite is destined for cell respiration and biomass (i.e., TCA cycle, alanine synthesis), whilst the rest is exported to the co-culture medium.Soon after pyruvate is incorporated into the myeloma cell's metabolic network, it is destined to support macromolecule synthesis, which we suspect is mainly destined for fatty acids.
Indeed, we can show that in-vitro co-cultured BMMSCs and myeloma cells consume less glucose than when in mono-culture (Fig 7).These observations suggest that these two cell types work as a metabolic community that promotes more efficient nutrient use within the bone marrow ( Fig 7).Experiments also revealed that co-cultured HS-5 cells showed increased glucose-derived pyruvate, while little change was observed in myeloma cells (Fig 7).This could indicate that pyruvate might be trafficked between cell types, potentially resulting in the enhanced efficiency of glucose metabolism, suggesting the possible validation of the model's prediction.
The model's flux distribution exhibits that the ratio of pyruvate over lactate in HS-5 cells is asymmetric towards pyruvate, in agreement with our co-culture experiments.It is unclear, however, what is the destiny of the lactate.This is because the export and import of this metabolite are balanced with its synthesis to pyruvate and back, presumably to maintain the cytosolic redox balance of NAD + /NADH.
The model predicts that once pyruvate is inside the mitochondria, it might be catabolised mainly via pyruvate dehydrogenases (PDH) and to a lesser extent via the pyruvate carboxylases (PC), where it could follow the TCA cycle terminating as L-Citrate, catalysed by the citrate synthase (CS).At this point, the model indicates that the metabolite is extruded from the mitochondria via citrate transport proteins in exchange for α-ketoglutarate, which we suspect is destined to synthesise fatty acids.As this flux is relatively large, the model bypasses the conversion of citrate to isocitrate by the aconitate hydratases, resulting in a relatively low proportion of mitochondrial glucose-derived a-ketoglutarate.These observations agree with our experimental data, in which we saw a decrease in the m+2 isotopomer from citrate to a-ketoglutarate (Fig 5B ), indicating that mitochondrial a-ketoglutarate in BMMSCs cells, in its majority, is not glucose-derived.The model predicts that this extrusion aids the balance of αketoglutarate fluxes from mitochondrial citrate exchange at a steady state -previous studies of inhibition of aconitate hydratases due to high levels of ROS saw similar results [95].According to our simulations, the BMMSC's TCA cycle, beyond the L-citrate step, is rescued via anaplerotic reactions due to high glutaminolytic activity, a hallmark of cancer (Fig 5D).Finally, in our in-silico model, we also observed a considerable flux by the L-lactate dehydrogenases (LDH) from glucose-derived pyruvate, indicating that the BMMSC is also highly glycolytic; as expected (The Warburg effect).This glucose-derived L-lactate is then extruded.However, we cannot rule out cells re-import this lactate once extruded.The result may support previous studies where these cells observed a paracrine-like behaviour, where L-lactate is reabsorbed [28].
Following the predicted incorporation of pyruvate into the mitochondria by the myeloma cell, our in-silico model suggests high rates of L-Citrate production, a possibility corroborated by our experimental results ( Fig 5).As in the BMMSC, pyruvate-derived L-Citrate is mainly exported to the cytosolic compartment in exchange for a-ketoglutarate (specifically via the mitochondrial carrier CTP-Slc25a1).Once in the cytosol, L-Citrate is mainly converted into fatty acids; the model's flux distribution selects long-chain fatty acids (i.e., palmitate).These fatty acids are primarily destined for biomass production, while a relatively small proportion is recycled back into the mitochondria via β-oxidation and incorporated back into the TCA cycle (Fig 6).Our model predicts that as most glucose-derived L-Citrate might be exported, this cell could exhibit low fluxes of isocitrate present in the TCA cycle of the myeloma model.Remarkably, these are even lower than those observed in the BMMSC cell model.Within the myeloma cell GEM, most of the available pyruvate is destined for the mitochondria and transported In co-culture, BMMSC produces relatively large amounts of pyruvate, with a fraction destined for cell respiration and the rest exported to the co-culture medium.Soon after, this pyruvate is uptaken by myeloma cells and incorporated into its carbon metabolism axis, destined to support biomass synthesis.However, both cells display a fractionated TCA cycle, in which L-citrate is exported, primarily for fatty acid synthesis, and the cycle is rescued via glutaminolysis.The model exhibits significant acetylation of asparagine to form N-acetyl-l-asparagine in myeloma cells.The metabolite is exported and uptaken by bone marrow mesenchymal stem cells, where reverse acetylation produces asparagine and acetate.This is used by the bone marrow mesenchymal stem cell for exogenous serine exchange, which is pipelined, in part, to produce pyruvate, lipids, antioxidants, and biomass.https://doi.org/10.1371/journal.pcbi.1011374.g007 passively by the mitochondrial membrane pyruvate carriers (MPCs).However, an intriguing result of the model's flux distribution describes that a significant proportion is transaminated into L-alanine in the cytosol and then transported into the mitochondria via an unspecified mitochondrial transporter.Once inside the mitochondria, L-alanine is de-transaminated back into pyruvate.The result is intriguing, as L-alanine is also a biomass precursor.Similarly to the BMMSC, most of the α-ketoglutarate necessary to resume the TCA cycle is supplied via the conversion of mitochondrial glutamate to α-ketoglutarate by L-Glutamate dehydrogenases.Finally, unlike the BMMSC results, the myeloma cell's flux distribution predicts that besides obtaining α-ketoglutarate from L-Citrate transport and glutaminolysis, the cell balances its mitochondrial α-ketoglutarate requirements by increasing the flux through the mitochondrial α-ketoglutarate/malate transporter (expressed by the Slc25a10/11 genes).This serves a double purpose, contributing further to sustaining anabolism and the redox balance needed to maintain the cell's high glycolytic demand.
Finally, we investigated additional possible intercellular metabolic exchange nodes across the multicellular model.Our model suggests the possibility of significant acetylation of asparagine to form N-acetyl-l-asparagine in myeloma cells (Fig 6).The metabolite is exported and uptaken by the BMMSC model, where reverse acetylation produces asparagine and acetate.In our model, this is used to exchange it for extracellular serine, which is destined, in part, to produce pyruvate, lipids, antioxidants, and biomass.As the myeloma cells show elevated levels of glutamine production for biomass synthesis and energy requirements, they might not be able to proliferate without exogenous asparagine.This indicates a demand which exceeds what can be generated through the GS/ASNS pathway.Asparagine, therefore, is N-acetylated using acetyl-CoA as the acetyl-donor to form N-acetyl-l-asparagine.In this way, the cell uses N-acetyll-asparagine derived-asparagine as a factor to regulate biomass synthesis and, perhaps, cell proliferation.The presence of this behaviour in our model is remarkable, as asparagine has been theorised to be a critical metabolite for proliferating cancer cells.

Discussion
Tumour metabolism remains an exciting space in which to identify novel drug targets.However, studies of cancer cell metabolism alone run the significant risk of identifying redundant metabolic enzymes in the cancer cell due to a more comprehensive metabolic network within tumours, including the neighbouring cells [96].Understanding the physiology of these metabolic features promises the creation of a targeted approach wherby intercellular metabolic cross-talk also represents an exciting space for therapeutic manipulation, restricting the tumoural metabolic network and decreasing the resilience of the cancer cell.
Although reprogrammed metabolic networks are not fully understood due to a myriad of experimental limitations, we can provide insights into the cellular metabolic status of cells through in-silico models [97][98][99].The ability of data integration, along with cell specificity in these models, may provide a broader picture, proving an invaluable tool to experimentalists.Furthermore, our proposed pipeline (Fig 2) is also relevant in studies where modelling cell metabolism within the context of the tumour microenvironment, which constitutes a complex adaptive ecological system, is critical, and due to the difficulty in experimentally probing the interactions between multiple cell types remains challenging [99].Such methodology (Fig 2 ) presents an alternative to known algorithms that sacrifice data integration for reductionist techniques.
Our proposed pipeline provides predictions addressing a known gap in the field of multiple myeloma.We reconstructed a state-of-the-art co-culture in-silico model to probe the tumour metabolism at the genome-scale whilst providing the necessary granularity to underpin the intercellular metabolic cross-talk between cancer and stroma in the context of the framework for constraint-based modelling.In our simulations, we predict that the maximal in-silico growth rates in both mono-culture and co-culture were consistent with the experimental counterparts (Fig 3A and 3B) and the findings reported in the literature [65].However, there is a nuance that merits attention: the growth rate of the HS-5 model closely aligns with the MM experimental data, and, interestingly, the MM model's growth rate mirrors that of the HS-5 experimental data.At first glance, this might seem inconsistent, yet it's crucial to remember that both models' growth rates fall within the range of experimental error.This detail is significant, as the experimental data does not categorically differentiate these cell types based solely on their growth rate, suggesting that growth rate may not serve as a robust, standalone validation criterion for cell-type specific modelling.Furthermore, the potential influence of numerous factors, including the complexity of in-vivo conditions and the culture medium will impact these rates.Assuming a maximal growth rate might be reductionist, given that a cell population's objectives might shift in response to the environment [14].Consequently, whilst we do not treat growth rate as the only validating metric for our model, it is an integral component that aligns with our observations.As a consequence we also recognise the importance of expanding our experimental data to include various distinct mediums to more accurately assess our model's predictions, given the current lack of regulatory constraints in our model [14,100].Thus, in future studies, we plan to incorporate proteomic and metabolimic data to help validate our model's predictions, which has demonstrated considerable potential for enhancing the support for experimentally observed phenomena, particularly when modelling the effects of perturbations in cell metabolism [100].
Taking the model's predictions into account, we carried out what we refer to as "phenotype tuning".This consists in integrating 13 C labelled glucose and glutamine data to constrain the model and produce physiologically relevant constraints.Our model predicts that myeloma cells may exhibit an effect characterised by the malignancy's potential reliance on glycolysis, leading to significant amounts of lactate production-this would be consistent with many previous studies on multiple cancer cell types [101,102].For instance, it has been hypothesised that lactate is not only a by-product of a malignancy's highly glycolytic metabolism but can be salvaged by other cell types for mitochondrial oxidative energy production [28,103,104].Similarly, it has also previously been suggested that lactate could be supplied by nearby stromal cells for the same reason-the directionality depending on the relative redox state of each cell [28,104].Furthermore, our model predicts that when in co-culture, BMMSCs and myeloma cells (HS-5 and JJN-3, respectively) form a metabolic community spanning across the tumour microenvironment where instead of lactate, as previously thought, pyruvate is trafficked from stroma to myeloma and not in the other direction (Figs 6 and 7).This intercellular interaction may be arbitrated unequivocally by transporters with a higher affinity for pyruvate over lactate, thereby establishing that these mechanisms rely mainly on pyruvate and not lactate as previously hypothesised [28,103].Interestingly, our computational model predicts that the pyruvate imported by myeloma cells could potentially support the cell's anabolism and energy generation.If our in-silico model predictions hold up under experimental testing, the result suggests that an alternative mode of action may be at play, which is aimed at supporting cell viability under limited resource conditions, which leads to an optimised metabolic community.Furthermore, with pyruvate transport targeting, lactate kinetics could be a powerful modality to perturb and increase the myeloma cell's vulnerability to conventional drug therapies.Such a target comes timely, as historically, it has been challenging to find a suitable therapeutic window when targeting multiple myeloma metabolism because most agents focus on phenotypes common to other normal tissues.This is evident from the years of unsuccessful trialling glycolytic inhibitors, and the side effects observed with the current anti-metabolite therapies such as 5-fluorouracil [105,106].These agents suffer from a lack of tumour specificity, often resulting in significant side effects [105,106].
A renaissance of research into cancer metabolism over the past two decades has been catalysed by significant improvements in our ability to detect and quantify thousands of different metabolites, alongside a revolution in sequencing technology that led to the identification of tumour-driving mutations in metabolic enzymes [107,108].Our study is part of a renewed drive towards developing novel agents targeting tumour metabolism.In this study we proposed an alternative in-silico reconstruction pipeline and used it to generate the first testable, integrated in-vitro/in-silico model of the metabolic network formed by malignant plasma cells and BMMSCs.To our knowledge, the majority of mathematical models that have been developed based on in-vitro data were produced under experimentally controlled conditions, where the metabolism of the malignant plasma cell in MM has characterised cell lines growing in isolation.Despite this, we are more likely to define efficacious cancer-specific targets for novel therapy development if the experimental system better re-capitulates the environment of the BM.Our integrated multidisciplinary workflows that iteratively utilise bench-side research alongside mathematical modelling can be utilised to identify and test the most appropriate targets for future pharmaceutical intervention.Introduction of metabolomics, fluxomics and growth-related data along with other related omics data into the reconstruction and validation process of GEMs can help increase the accuracy and prediction power of any model generated using our pipeline [109].Our workflow, as demonstrated, can be used to generate useful in-silico models that underpin difficult-to-see behaviours and predict their results alongside in-vitro experiments.

Fig 3 .
Fig 3.Growth rate comparison between models and experimental results in mono-culture and co-culture. A. Depiction of a comparative analysis of the growth rates of in-vitro and in-silico models.The figure illustrates the experimental growth of BMMSCs in-vitro culture (depicted as black dots), normalised for comparison.This growth curve is contrasted with the projected growth as determined by the exponential growth model that uses the growth rate predicted from the BMMSC Constraint-Based Model (CBM); solved utilising the Flux Balance Analysis (FBA) procedure with the maximum biomass production rate serving as the optimization objective.This is represented by the black intermittent line.Simultaneously, the figure demonstrates the measured, normalised growth of the MM cell line in-vitro culture, symbolized by red dots.This is set against the growth curve predicted, with maximum biomass production rate being the primary optimisation goal; denoted by the intermittent red line.B. Illustration of the growth rate comparison between the model and experimental results in co-culture.The figure demonstrates the measured, normalised in-vitro co-culture growth of the BMMSC cell line; (black dots).This is contrasted with the growth curve given by the exponential growth model solution, which is fitted with the growth rate predicted by the in-silico co-culture model, represented by the intermittent black line.Concurrently, the figure presents the measured, normalised in-vitro co-culture growth of the MM cell line, depicted by red dots.This is set against the growth curve yielded by the exponential growth model solution, adjusted to the growth rate forecasted by the in-silico co-culture model, denoted by the intermittent red line.Note that in both cases, the monoculture and the co-culture, the negligible differences between the modelled and measured growth rates of BMMSCs and MMs limit cell type discrimination.

Fig 4 .
Fig 4. Mitochondrial respiration experiment and in-silico simulation using the Oroboros respirometer to measure mitochondrial respiration capacity.Panel A. contrasts the results of an Oroboros respiration study carried out on bone marrow mesenchymal stem cells (BMMSCs) with predictions from our in-silico model.The superimposed illustration provides a clear comparison, demonstrating the extent of congruence between the model's predictions (red) and the experimentally observed data (purple).Panel B. sets the results of the Oroboros respiration study on the myeloma cell line against the predictions of our in-silico model.The overlaid results allow for an explicit comparison, showcasing the correlation between the model's predicted outcomes (yellow) and the data derived from experimental observation (green).https://doi.org/10.1371/journal.pcbi.1011374.g004 Fig 5 summarises the results of our metabolic flux analysis.Panels A and C of Fig 5,

Fig 5 .
Fig 5. 13 C-MFA simulations.A. & B. Mass isotopomer distributions from [U-13 C]Glucose and [U-13 C]Glutamine tracing.The Mass Isotopomer Distributions represent the relative abundance of each isotope, normalised to the total sum of all possible isotopologues.Glucose-derived metabolite abundances.Data suggests high glycolytic activity in both cells.In JJN-3 cells, pyruvate-derived alanine is much higher than in HS-5 cells.The abundance of pyruvate in JJN-3 cells is also more significant than that of HS-5.It is possible to see that although pyruvate-derived citrate follows the expected TCA cycle, it does not continue downstream as aketoglutarate is low in JJN-3 cells.This indicates an interruption of the TCA cycle.C. & D. Mass isotopomer distributions from [U-13 C]Glucose and [U-13 C] Glutamine tracing.Glutamine-derived metabolite abundances.Data suggests that the TCA cycle in both cells is supported via anaplerosis, where exogenous glutamine is interconverted to a-ketoglutarate and then incorporated into the TCA cycle.

Fig 6 .
Fig 6. A. Incorporation of glucose carbons into pyruvate is also altered by co-culture; while HS-5 cells show lower pyruvate incorporation in mono-culture compared to co-culture conditions, there is little change in JJN-3 incorporation.B. When incubated with 13 C 6 Glucose, label enrichment into the TCA cycle metabolite malate in both the myeloma cell line (JJN-3) or BMMSCs (HS-5) indicates a significant increase in oxidative glucose metabolism (m +2 isotopomer) as well as anaplerotic use (m+3 isotopomer).C. Co-culture of BMMSCs and MM cell lines, HS-5 and JJN-3, respectively, results in reduced glucose use compared to either cell type in mono-culture (corrected to cell number), suggesting a more efficient use of glucose carbons when both cells are cultured together.(In-figure referencing: ***; p<0.001, ****;p<0.0001.From Anova and Dunn's multiple comparison).https://doi.org/10.1371/journal.pcbi.1011374.g006

Fig 7 .
Fig 7. Model predictions.In co-culture, BMMSC produces relatively large amounts of pyruvate, with a fraction destined for cell respiration and the rest exported to the co-culture medium.Soon after, this pyruvate is uptaken by myeloma cells and incorporated into its carbon metabolism axis, destined to support biomass synthesis.However, both cells display a fractionated TCA cycle, in which L-citrate is exported, primarily for fatty acid synthesis, and the cycle is rescued via glutaminolysis.The model exhibits significant acetylation of asparagine to form N-acetyl-l-asparagine in myeloma cells.The metabolite is exported and uptaken by bone marrow mesenchymal stem cells, where reverse acetylation produces asparagine and acetate.This is used by the bone marrow mesenchymal stem cell for exogenous serine exchange, which is pipelined, in part, to produce pyruvate, lipids, antioxidants, and biomass.